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ABSTRACT 

We describe the implementation of the shearing box approximation for the study of the dynamics of accretion 
disks in the Athena magnetohydrodynamics (MHD) code. Second-order Crank-Nicholson time differencing is 
used for the Coriolis and tidal gravity source terms that appear in the momentum equation for accuracy and 
stability. We show this approach conserves energy for epicyclic oscillations in hydrodynamic flows to round-off 
error In the energy equation, the tidal gravity source terms are differenced as the gradient of an effective potential 
in a way which guarantees that total energy (including the gravitational potential energy) is also conserved to 
round-off error. We introduce an orbital advection algorithm for MHD based on constrained transport to preserve 
the divergence-free constraint on the magnetic field. This algorithm removes the orbital velocity from the time 
step constraint, and makes the truncation error more uniform in radial position. Modifications to the shearing box 
boundary conditions applied at the radial boundaries are necessary to conserve the total vertical magnetic flux. In 
principle similar corrections are also required to conserve mass, momentum and energy, however in practice we 
find the orbital advection method conserves these quantities to better than 0.03% over hundreds of orbits. The 
algorithms have been applied to studies of the nonlinear regime of the MRJ in very wide (up to 32 scale heights) 
horizontal domains. 

Subject headings: hydrodynamics, MHD, methods:numerical 



1. INTRODUCTION 

' Numerical methods are an important tool for studying the 
nonlinear gas dynamics in accretion flows. For example, much 
of what we have learned about the saturation of the magne- 
, torotational instability (MRI) has come from magnetohydrody- 
namic (MHD) simulations (Balbus & Hawley 2003). A large 
, fraction of such calculations adopt a local approximation (Hill 
' 1878) termed the "shearing box", first introduced in studies of 
■ the MRI by Hawley et al. (1995, hereafter HGB). In this ap- 
proximation, the equations of motion are written in a local, 
Cartesian reference frame co-rotating with the disk at some ar- 
. bitrary radius ro. The approximation is valid provided the linear 
extent of the domain under study is small compared to ro. The 
shearing box approximation has limitations (Regev & Umurhan 
2008), in particular it cannot be used to calculate important 
global properties of the disk like the net mass accretion rate, ra- 
dial density and temperature profiles, or the spectrum of emitted 
radiation. Nevertheless it has provided a useful laboratory for 
the study of important questions related to the local dynamics 
of accretion disks (Balbus 2003). For example, recent studies 
that use the shearing box include the effect of microscopic dif- 
fusivities on the saturation level of MHD turbulence driven by 
the MRI (Fromang et al. 2007; Lesur & Longaretti 2007; Si- 
mon & Hawley 2009), the saturation of the MRJ in radiation 
pressure dominated disks (Hirose et al. 2008), and the effects 
of non-ideal MHD and dust on the properties of the MRJ in pro- 
tostellar disks (Turner & Sano 2008; Ilgner & Nelson 2008). 



The shearing box approximation requires source terms be 
added to the momentum and energy equations, and special 
boundary conditions be used in the radial direction. Most of the 
simulations of the MRJ presented to date have used operator- 
split methods like ZEUS (e.g. HGB; Stone et al. 1996), or 
simple finite-difference methods like the PENCIL code (e.g. 
Johansen et al. 2009). In both cases, adding the shearing box 
source terms is straightforward (HGB). More recently, higher- 
order Godunov methods that adopt the conservative form have 
begun to be applied to studies of accretion flows in the shearing 
box (Shen et al. 2006; Fromang & Papaloizou 2007; Bodo et 
al. 2008; Piontek et al. 2009; Simon et al. 2009; Tifley et al. 
2009; Gressel 2010). 

In this paper, we provide a detailed description of the al- 
gorithmic extensions for the shearing box approximation in 
Athena, a recently developed higher-order Godunov code for 
astrophysical MHD. The basic MHD algorithms in Athena are 
documented in Gardiner & Stone (2005a; 2008), and details of 
the implementation and tests of the methods are given in Stone 
et al. (2008, hereafter SGTHS) and Stone & Gardiner (2009). 
Currently there are two versions of Athena, one implemented 
in C and the other in Fortran. An extension of the Fortran ver- 
sion for the shearing box has already been used for new stud- 
ies of the MRI (Simon & Hawley 2009; Simon et al. 2009). 
The numerical algorithms in the C version are slightly differ- 
ent. They were first introduced by Gardiner & Stone (2005b), 
and used to study hydrodynamic shearing waves by Shen et al. 



2 



(2006) and Balbus & Hawley (2006). They have the advan- 
tage of preserving the energy integral in epicyclic motion to 
round-off, conserving the total (including gravitational poten- 
tial) energy to round-off, and producing virtually no aliasing of 
trailing into leading waves. We emphasize, however, that since 
epicyclic motion is destabilized by a weak magnetic field (Bal- 
bus & Hawley 1991), the accuracy of this approach can only 
be quantitatively demonstrated for hydrodynamic flows. Al- 
though similar methods were adopted and extended by Gressel 
& Ziegler (2007), this paper provides the first comprehensive 
description of the algorithms in the C version of Athena. The C 
version is also being used for new studies of the MRI (Davis et 
al. 2010). 

In addition to methods for the shearing box source terms, we 
also describe the implementation of an orbital advection algo- 
rithm for MHD in Athena. Orbital advection can greatly in- 
crease the efficiency of calculations in domains that span more 
than one scale height in the radial direction (so that the differ- 
ence in the orbital speed across the domain is supersonic), since 
the time step constraint for stability does not depend on the 
magnitude of the local orbital velocity, but only on the ampli- 
tude of the fluctuations in the velocity around this value. More- 
over, orbital advection can improve the accuracy of the inte- 
gration by making the truncation error more uniform in radius 
(Johnson et al. 2008, hereafter JGG). Orbital advection was first 
introduced for hydrodynamic studies of disks by Masset (2000) 
in the FARGO code. More recently, JGG have described an ex- 
tension of orbital advection methods to MHD for ZEUS -type 
codes. Johansen et al. (2009) have described the extension of 
the PENCIL code with orbital advection using Fourier trans- 
form methods. The method we have implemented in Athena 
is quite different from these previous approaches. In particu- 
lar, we update the magnetic field using constrained transport 
(Evans & Hawley 1988) to guarantee the divergence-free con- 
straint is enforced to machine precision, using an effective elec- 
tric field produced by the orbital motion. This greatly simplifies 
the method. 

The organization of this paper is as follows. In the following 
section, we catalog the basic equations solved by Athena in the 
shearing box approximation. In §3 we describe our implemen- 
tation of the source terms in these equations, including tests of 
our methods. In §4 we describe the shearing box boundary con- 
ditions for the conservative variables, including issues associ- 
ated with parallelization of the algorithms with MPl. Finally, in 
§5 we describe our orbital advection method, along with tests. 
We present preliminary simulations of the MRI in very wide 
horizontal domains and summarize in §6. 
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where T is the total stress tensor 

T = (/:'+bV2)I-BB, (5) 

p is the gas pressure, E is the sum of the internal, kinetic, and 
magnetic energy densities 
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and = V • V, = B • B. The shear parameter q is defined as 

1 d\nn'^ 



q = - 



(7) 



2 d\m 

so that for Keplerian flow ^ = 3/2. The other symbols have their 
usual meaning. These equations are written in units such that 
the magnetic permeability p= I, and an equation of state ap- 
propriate to an ideal gas has been assumed in writing equation 
121 that is = (7- l)e (where 7 is the ratio of specific heats, and 
e is the internal energy density). The techniques described in 
this paper are easily generalized to other equations of state. 

The first two source terms on the RHS of equation|2]represent 
the radial and vertical components of the gravitational force in 
the rotating frame, while the third term represents the Coriolis 
force. In the energy equation [3] the two source terms on the 
RHS represent the work done by the radial and vertical compo- 
nents of gravity in the rotating frame. Most of the challenges 
associated with the addition of the shearing box source terms 
to Godunov schemes like Athena are related to the tidal gravity 
and Coriolis terms which act in the orbital {x-y) plane, there- 
fore in the discussion below we describe the algorithms for the 
source terms in the horizontal and vertical directions separately. 

The equations of motion in the shearing box admit a simple 
equilibrium solution representing uniform orbital motion, 

\K = -q^()xS- (8) 

Note this velocity is time-independent, varies only in the 
x-direction, and only the y-component is non-zero. These 
properties are important for the orbital advection scheme de- 
scribed in §5. 



2. BASIC EQUATIONS 

The local shearing box approximation adopts a frame of ref- 
erence located at radius ro corotating with the disk at orbital 
frequency f7o = f^('"o)- In this frame, the equations of MHD are 
written in a Cartesian coordinate system (x,y,z) that has unit 



3. SHEARING BOX SOURCE TERMS 

One approach to implementing source terms in Godunov 
schemes like Athena is to operator split them from the flux di- 
vergence terms, resulting in a system of ODEs that can be inte- 
grated with any number of methods. For example, Simon et al. 
(2009) have successfully used this approach to study the satu- 
ration of the MRI with the Fortran version of Athena. However, 
in studies of decaying hydrodynamic turbulence in the shearing 
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box (Shen et al. 2006), we found the kinetic energy in velocity 
fluctuations could increase at very late times, rather than decay 
monotonically to zero. We traced the source of this growth to 
an artificial amplification of epicyclic motions by the trunca- 
tion error associated with the source terms in the momentum 
equations. As we show below, in hydrodynamics it is possi- 
ble to develop a discretization that will instead conserve this 
energy to round-off error While this approach may not pro- 
vide significant improvement for MHD flows, since epicyclic 
motion is unstable with weak magnetic fields, we show below 
it is more accurate for hydrodynamics, and therefore we have 
adopted it as the default algorithm in both hydrodynamics and 
MHD. The methods for integrating the source terms described 
below were first introduced by Gardiner & Stone (2005b), sub- 
sequently Gressel & Ziegler (2007) have adopted and expanded 
upon this approach. 

To describe our methods, it is useful to define the x- and 
y-components of the momentum density fluctuations 

m = pv^ (9) 

11 = p(Vy + qilox) (10) 

To investigate the properties of epicyclic motion in a stress-free 
medium, we consider the case in which the velocity fluctuation 
is a function of time and rewrite the momentum equation|2]us- 
ing these variables, giving 



dm 

— = 2Qon 
ot 

dn 

— = (q-2)nom 
ot 



(11) 
(12) 



Multiplying equationfTTIbv m, equation [T2]bv -2n/{q-2), and 
adding gives 

5 / 1 2 



dt 



+ ■ 



2 — q 



= 



(13) 



Thus, there is a conserved energy which for Keplerian flow 
iq = 3/2) is £epi = [m2 + 4n2]/(2p). 

We would like the numerical discretization of the momen- 
tum equation to conserve Egpi exactly. Remarkably, nothing 
more complicated than Crank-Nicholson time differencing is 
required. Following the usual convention, we discretize time 
into non-uniform steps, and use a superscript to denote the 
time level of any quantity, with the timestep St defined as 
St = t"^^ -t". We define time difference and averaging opera- 
tors as 

(14) 



x=^-{x"^'+x") 



(15) 



respectively. Then, the Crank-Nicholson time difference for- 
mulae for equations [TTI and [T2l are 

[m]=2noStn (16) 
[n] = (q-2)noStm (17) 

By multiplying the first of these equations by m, the second by 
-2n/(q-2), and adding gives, after some manipulation 

[^2^J_„2] = 0. (18) 
1 — q 



Thus, Crank-Nicholson discretization conserves the energy in- 
tegral for epicyclic oscillations in a discrete sense. There is a 
simple physical interpretation of this result. The Coriolis force 
Fc = (-2f2 X v), hence v • Fc = and the Coriolis force can do no 
work. Without care, the discretized Coriolis force may not be 
orthogonal to the velocity (in a time averaged sense). This can 
lead to unphysical growth or decay of the energy in epicyclic 
motion. Our tests have revealed that a forward Euler discretiza- 
tion leads to a growing amplitude for epicyclic oscillations, 
while a backward Euler discretization leads to a decaying am- 
plitude. Fortunately, the average of the two (Crank-Nicholson 
differencing), conserves the energy. Similarities can be drawn 
between integrating fluid motion due to the Coriolis force, and 
particle orbits in a central potential, e.g. forward Euler applied 
to the latter leads to outward spiraling rather than closed orbits. 

3.1. Source Terms in the Momentum Equation 

We now develop a finite-volume discretization of the mo- 
mentum equation |2] using a Crank-Nicholson time discretiza- 
tion of the tidal gravity and Coriolis source terms. We con- 
sider only the x- and y-components of this equation. Written 
in terms of the momentum density fluctuations, m and n respec- 
tively, these equations are 

0171 

■ + V-F,„ = 2r!on (19) 



dn 
'dt 



dt 

+ V-F„ = (^-2)f^o'« 



(20) 



where F^ and F„ are vectors whose components are the fluxes 
of the momentum density fluctuations in each direction. These 
are related to the vector of fluxes of the corresponding com- 
ponents of the momentum density, Fp,,^ and Fp,,, respectively, 
via 

F,„ = Fp„, (21) 

F„ = Fp,,, + ^l]axFp (22) 

We now adopt a finite- volume discretization of equations [T9l 
and l20l (see §3 in SGTHS for a discussion of this approach). 
We use indices ii,j,k) to denote spatial locations on a discrete 
grid with cell centered locations {xi,yj,Zk)- Half integer indices 
are used to denote the appropriate cell interfaces. The depen- 
dent variables stored on this mesh m'jj/^ and n'jj/^ are under- 
stood to be volume averaged values in the sense of equation 12 
in SGTHS. Integrating over a timestep St = f""*"' - f" and over the 
volume of a cell, equations [T9]and|20lbecome 



n+l 



St 



+ (V-F„,),.^., = 2r!o 



(23) 



(mf^.^ + m^'y^) 

+ (V • F„),,. , = {q- 2)»o ^ (24) 

In a finite volume approach, the time- and volume-average of 
the flux divergence (the second term in both equations) would 
be rewritten, using the divergence theorem, as the difference 
of the time- and area-averaged fluxes at each of the cell faces, 
e.g. equations 11-15 in SGTHS. In a Godunov scheme, these 
fluxes are computed with a Riemann solver. For brevity, we 
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have suppressed expanding the flux divergence terms into these 
differences, with the understanding that our notation is meant 
to represent these terms. 

Solving equation l24l for and substituting into 

equation|23]gives, after some manipulation 



<)Ji = '«">..-'5f(V-F„,),^., 



2 + i2-q){floSty 



dt- 



-(V-F„),^,, 



(25) 



Similarly, solving equation |23] for (m" j i^ + mlyj^) and substitut- 
ing into equationl24l gives 

nij'k = <ik-StW^),jM 



(26) 



2(q-2)noSt \ (j „ -JfT^y-p^ \ 
+ (2-,)(llo^r)0 lr--'^"2^^-^"'^--^4 



In order to convert these update equations for m, j ^ and n, yjt 
respectively, into update equations for {pvx)i,j,k and {pv^)i jk 
(the discrete form of the conserved quantities actually updated 
in the code), relationships between the discrete fluxes of these 
quantities are required. For any cell (i,j,k) with zone-center 
jc-position x,-, the relationships between the components of the 
volume averaged momentum density fluctuations and the mo- 
mentum density are 

mjk = (pvx)ijk (27) 

nijk = {pvyjijk + q^apijkXi (28) 

From equations |2T| and |22] the divergence of the fluxes of these 
quantities are related via 

V-F„, = V-Fp., (29) 

V • F„ = V • Fp„,, + ^r^oV • (xFp) (30) 

which can be written in a finite volume discretization as 

W^)i.j,k = W^i,j.k (31) 



(V • F„),^. , = (V • Fp„) . .^^ + qnm (V • Fp) . .j^ 

+ {q/2)nm {iFp\i-xii + (^'p)x,+i/2) (32) 
Finally, inserting equations |28] and l32]into |26] gives 



{pvyyi:^^, = {pv.)lj,,-5t{V-¥^,,,,)..^ 



2(q-2)noSt 
2 + {2-q)(nodt)^ 



St 



"j.k~ '^m)ij,k 



(33) 



Since the m/./jt and (pvx)ij,k are identical, the update relation 
for the latter is given by equation |25] 



Equations |25] and |33] above represent the desired finite vol- 
ume discretization of the momentum equation, where the diver- 
gence of the fluxes of the momentum fluctuations that appear 
in these equations are given in terms of the divergence of the 
fluxes of the conserved quantities actually returned by the Rie- 
mann solver by equations [311 and [32l As we will show through 
the tests described in §3.4, the Crank-Nicholson time differ- 
encing of the shearing box source terms used in these equations 
conserves the energy integral associated with epicyclic motion 
exactly. 

To implement these difference equations in a computer code 
additional algorithmic steps are required. For example, the 
reconstruction of the left- and right-interface states in the 
jc-direction (see §4.2 of SGTHS) which is needed to compute 
the time-and area averaged fluxes of the conserved variables at 
the jc-interfaces using a Riemann solver requires the addition 
of shearing box source terms. In particular, before the left- and 
right-states in the primitive variables at x-interfaces are con- 
verted back to the conserved variables (q/^ ,_i/2 and qs ,_i/2 in 
the notation of SGTHS) at the end of the first-, second-, or third- 
order reconstruction algorithms describe in §4.2 of SGTHS, the 
shearing box source terms for the velocity components Vx and 
Vy must be added for 5t /2. We have found that simple forward 
Euler time-differencing is adequate for this step. 

Similarly, the directionally unsplit corner transport up- 
wind (CTU) integrator used in Athena (see §5.1 and §6.1 in 
SGTHS) uses transverse flux gradients to correct the inter- 
face states in multidimensions. For the flux gradients in the 
x-direction (added to the y-interface states in 2D, and the 
y- and z-interface states in 3D), the appropriate shearing box 
source terms must be added to x- and y-components of the mo- 
mentum for 5t/2. Again, we have found that simple forward 
Euler differencing is adequate for this step. 

Finally, in order to compute the cell-centered reference elec- 
tric field at the half time step needed for the CT algo- 
rithm (computed in step 5 of the 2D CTU algorithm, or step 6 
in the 3D CTU algorithm), the shearing box source terms with 
forward Euler discretization must be added in the calculation of 
the velocity components at the half time step. 

In summary, our algorithm for the momentum equation up- 
date in the shearing box approximation as implemented in 
Athena consists of the following modifications to the CTUh-CT 
algorithm described in SGTHS: 

1. Add shearing box source terms to the left- and right- 
states for the velocity at x-interfaces only during recon- 
struction step, before converting the reconstructed prim- 
itive variables to conserved variables. 

2. Add shearing box source terms to the x- and 
y-components of the momentum when a transverse flux 
gradient in x-direction is applied to y- and z-interface 
states as part of CTU integrator. 

3. Add shearing box source terms to velocity when com- 
puting the cell-centered reference electric field at the 
half time step needed by the CT algorithm (see step 5 
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in §4.2 of Gardiner & Stone 2008). 

4. Use equations|25]and[33]to evolve the momentum in the 
final conservative update, where the divergence of the 
fluxes of the momentum fluctuations come from equa- 
tion s|3l] and [32] 

As we show with tests, the main advantage of this algorithm 
is that it conserves the energy integral for epicyclic motion ex- 
actly. 

3.2. Source Terms in the Energy Equation 

Next we consider the finite-volume discretization of the en- 
ergy equation [3] including the shearing box source terms. The 
key to this discretization is to define an effective potential for 
the shearing box 

^SB = -q^V+\^V (34) 

so that the source terms in the energy equation[3]can be written 
as 

dE 

— + V- [£'v + T-v] = pv- V$sB (35) 
at 

When written in this form, it is clear that the source term repre- 
sents the rate of change in the gravitational potential energy per 
unit volume. 

We seek a discretization of this term such that when inte- 
grated over volume, we recover the rate of change of energy 
due to the work done at the boundaries. This suggests that for 
each computational cell, we discretize the source term by the 
difference of the work done at the edges of the cell. Ignoring for 
the moment the vertical component of gravity, so that the effec- 
tive potential depends only on x, the appropriate finite-volume 
discretization of equation[35]is 

dt ( ^SB,i+\ I2,i,k - ^SB,i,j,k 
+ (^'p).v,,-l/2,-,. ^ j (36) 

where (fp)A.(-i-i/2.>.ji: are the mass fluxes in the x-direction re- 
turned by the Riemann solver at the x-interfaces. It is straight- 
forward to confirm that when integrated over volume, the 
source term reduces to the net mass flux times the difference 
in the gravitational potential across the domain. Thus, the only 
route via which the total energy in the domain can change is 
through work done at the boundaries. Within the domain en- 
ergy can be exchanged between its kinetic, magnetic, and ther- 
mal forms, but it cannot be lost as truncation error Even though 
we evolve an energy variable that does not contain the gravita- 
tional potential energy, we use a discretization that conserves to 
machine precision the total energy, that is E + p^sB- 

In addition to the use of the update equation |36] one other 
algorithmic step is required to include the shearing box source 
terms in the energy equation. This step is to include source 
terms to the transverse flux gradient corrections to the left- and 



right-interface states for E in the multidimensional CTU al- 
gorithm. To be more precise, when the flux gradients in the 
x-direction are added to the y-interface states of E in 2D, 
and the y- and z-interface states of E in 3D, the appropriate 
shearing box source terms must included for 6t /2. Again, we 
have found that simple forward Euler differencing is adequate 
for this step. Note that since the reconstruction step uses the 
primitive {p) rather than conserved (E) variables, no additional 
source terms are required for the evolution of the energy equa- 
tion in the reconstruction step. 

Thus, the extensions to the algorithms for energy equation 
update in the shearing box approximation consists of the fol- 
lowing modifications: 

1. Add shearing box source terms to E when a trans- 
verse flux gradient in x-direction is applied to y- and 
z-interface states of E as part of the CTU integrator. 

2. Use equation [36] to evolve E in the final conservative 
update. 

As shown in Gardiner & Stone (2005b) in studies of the MRI 
in the shearing box with Athena, with this algorithm the rate 
of change of E is equal to the work done at the boundaries to 
machine precision. 

3.3. Including Vertical Gravity 

To simplify the discussion in the previous subsections, we 
considered only the two components of the momentum equa- 
tion in the orbital (x-y) plane, and the source term associated 
with the vertical component of gravity in the energy equation 
was ignored. We have shown that care is needed in the dis- 
cretization of the source terms in these equations in order to 
correctly capture the dynamics of epicyclic motion. However, 
for the discretization of the gravitational acceleration term in 
the vertical (z-) component of the momentum equation, and the 
work term associated with the vertical component of gravity in 
the energy equation, we have found less complicated methods 
are adequate. 

In Athena, the vertical gravity terms in both the vertical com- 
ponent of the momentum and the energy equation are added by 
differencing the effective potential equation |34] For the verti- 
cal component of the momentum equation, several algorithmic 
steps are required: 

1 . At the end of the reconstruction step, vertical accelera- 
tion for 6t/2 is added to the left- and right-states of v, at 
z-interfaces using the vertical gradient of the potential. 

2. When the transverse flux corrections in the z-direction 
are applied to the x- and y-interface states as part of the 
CTU unsplit integrator, vertical acceleration for 6t/2 is 
added to pv. using the gradient of the potential. 

3. When computing the cell-centered reference electric 
field at the half time step needed by the CT algorithm 
(see step 5 in §4.2 of Gardiner & Stone 2008), vertical 
acceleration for St /lis added to v. using the gradient of 
the potential. 
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4. In the final conservative update, the vertical force for 5t 
is added to pv, using the gradient of the potential and the 
density evaluated at the half time step, ip" j ii + Pl^j\)/2. 

For the energy equation, the final conservative update equa- 
tion[36]must be extended to include a term Se-, where 

C ^^(tJ?\ ^SB,i,j,k+U2-^SB,iJ.k 

SEr. - -y [(Fphj.jMi/2 

+ (t'p)z.i.j,k-i/2 J (37) 

This leads to the following algorithmic steps: 

1. Add shearing box source terms to E when a trans- 
verse flux gradient in z-direction is applied to x- and 
y-interface states of E as part of the CTU integrator. 

2. Use equation [36] including the term given by equation 
l37]on the RHS, to evolve E in the final conservative up- 
date. 

Inclusion of the source terms for vertical gravity, along with 
appropriate boundary conditions in the vertical direction, allow 
studies of the MRI in stratified disks (e.g. Stone et al. 1996) to 
be continued using Athena (Davis et al. 2010). 

3.4. Tests of the Source Tenns 

There are a variety of axisymmetric {d/dy = 0) solutions to 
the MHD equations in the shearing box that serve as useful tests 
of the implementation of the source terms. Non-axisymmetric 
solutions required a more sophisticated treatment of the bound- 
ary conditions, as described in the next section. All of the tests 
in this section use a periodic domain in the (x-z) plane which 
spans -Lx/2 <x< L^/l and -L,/2 <z< The orbital fre- 
quency ilo = 10"^, the shear parameter ^ = 3/2, and an isother- 
mal equation of state is used with sound speed C, = 10"^. Un- 
less otherwise stated, the initial conditions consist of a uniform 
density medium with po = 1, pressure p = p^Cj, and orbital ve- 
locity Vy{x, 0) = -qilox. Third-order reconstruction, the HLLC 
(for hydrodynamics) or HLLD (for MHD) Riemann solver, and 
the CTUh-CT unsplit integrator are used in all the tests. 

A good first test is the evolution of epicyclic oscillations. We 
choose the initial radial velocity = O.IC,, where Cs = 10"^, 
in domains of size Lx = L~ = 1, 10, and 50, and evolve the flow 
for thousands of orbits using a grid of 64-. We vary the size 
of the domain while keeping the resolution constant in order 
to study the effect of very large timesteps on the accuracy of 
the integration algorithm. Figure 1 shows the oscillations in Vi 
over the first 20 orbits in each run. In each case, we observe 
epicyclic motion with constant amplitude: the energy integral 
in epicyclic motion is conserved exactly. Note there is a small 
dispersion error in the largest domain. However, in this case the 
time step is so large that there are only 14 timesteps per orbit. 
Thus, the small dispersion error (< 2%) observed in this case 
(due to the Crank-Nicholson differencing) is to be expected. In- 
creasing the resolution, and therefore decreasing the timestep, 
makes this error converge at second order. The key result is that 



even with large timesteps, the amplitude of oscillations is con- 
stant. However, we emphasize that this result is relevant only 
for hydrodynamics, since with weak magnetic fields epicyclic 
motion is unstable. 

A more quantitative test is provided by the propagation of 
nonlinear axisymmetric density waves in the shearing box (Fro- 
mang & Papaloizou 2007). The profile of such waves is given 
by the solution to two ordinary differential equations; for this 
test we use a numerical solution (kindly provided by S. Fro- 
mang) on a grid of 40 points in a domain with = 9.593 with 
a wave with amplitude Pmax/po = 1-05. Figure 2 shows the er- 
ror in the density, defined as 5p(x,t) = (p{x,t)- p{x,0))/ p(x,0) 
after f = 10, 50, and 100 wave periods T, where T = 5254.31. 
Note the error is smooth (there are no oscillations indicative of 
dispersion error), and never larger than 1.5%. We also plot the 
time evolution of the error in the wave amplitude, defined as 
(e- eo)/eo, where e = max{p)- po. After 100 wave periods, we 
find the amplitude has decreased by only 3%, which is a mea- 
sure of the low numerical diffusion rate in Athena. These re- 
sults compare very favorably to the test results presented in Fro- 
mang & Papaloizou (2007), which were computed at a much 
higher resolution. 

Finally, another useful test is the nonlinear evolution of ax- 
isymmetric modes of the MRI with no net flux. We use a 
domain of size = Lj = 1, and an initially vertical field with 
B = (0,0,i2po/ iSy^^sinlnx/Lx), where /3 = 4000. This test is 
identical to run Sic in Hawley & Balbus (1992). In figure 3, we 
show images of the azimuthal velocity fluctuations Vy + qilox at 
various times in the evolution. The growth and saturation of 
the MRI is evident, with transition to two-dimensional MHD 
turbulence at late times. In two dimensions, turbulence can- 
not be sustained by the MRI, and the magnetic energy even- 
tually decays. The rate of decay depends on numerical diffu- 
sion, which in turn depends on the grid resolution and numer- 
ical algorithm. Previously, we have used the rate of decay of 
the poloidal magnetic energy after saturation to compare the 
numerical dissipation rates of ZEUS and Athena (Stone & Gar- 
diner 2005; Stone 2009), to show that Athena gives lower de- 
cay rates at the same grid resolution. Figure 4 compares the 
time evolution of the energy associated with the radial compo- 
nent of the magnetic field B^./2 for various resolutions and both 
second- and third-order reconstruction. At every resolution the 
third-order reconstruction is measurably more accurate, in the 
sense that the rate of decay of the magnetic energy is decreased. 
Perhaps the best indicator of the accuracy of the methods is the 
evolution at very early times, while the amplitude of the MRI 
is very small. The short period of decay (due to the expansion 
of regions with the highest magnetic pressure), followed by ex- 
ponential growth over 10 orders of magnitude in energy from 
the initial (very small) amplitude indicates Athena has a very 
good dispersion relation for slow MHD waves (which was also 
shown in the Stone et al. 2008 linear wave convergence tests). 

4. SHEARING BOX BOUNDARY CONDITIONS 

As first discussed in HGB, non-axisymmetric solutions in the 
shearing box require special boundary conditions that offset the 
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solutions by the distance the radial edges of the domain have 
been displaced by the background shear Mathematically, the 
shearing box boundary conditions can be expressed as (HGB; 
Gressel & Ziegler 2007): 

f(x,y,z)^f{x±L,,yTwt,z) (38) 

where / = (p, pv^, pv,, B), and 

pvy{x,y,z)^ pVy{x±L^,y^wt,z)T (39) 

E{x,y,z) ^ E(x±L,,yTwt,z) T pVyW + pw^/l (40) 

where w = qftoLf is the difference in the orbital velocity across 
the domain. Although a variety of authors have already de- 
scribed the implementation of the shearing box boundary con- 
ditions in various codes (e.g. HGB; Gressel & Ziegler 2007), 
the details of the implementation can be important, thus we de- 
scribe our methods below. 

In Athena, we implement the shearing box boundary condi- 
tions by first applying periodic boundary conditions in the ra- 
dial direction, (with the appropriate shift in the y-momentum 
and energy given by equations [39] and |40] |. and then use a 
conservative remap of all quantities in the ghost cells in the 
y-direction by a displacement wt, the distance the boundaries 
have sheared in time t. This distance in general is not an integer 
number of cells. To remap the solution a fractional portion of 
a cell we use the same higher order reconstruction functions as 
used in the integrator to compute the amount of each conserved 
variable to assign to the remapped cells. For the magnetic field, 
we remap each component in the ghost zones directly, rather 
than evolving the field in the ghost zones using a remap of the 
emfs (e.g. Gressel & Ziegler 2007). This procedure does not 
guarantee the divergence free constraint is maintained in the 
ghost zones. However, since these cells are never used for any- 
thing more than interpolation, we have not found this to be a 
problem. 

As pointed out by Gressel & Ziegler (2007), the shearing box 
boundary conditions can destroy conservation, because the inte- 
gral of the fluxes of the conserved quantities over the two radial 
faces may not be identical due to the remap. In particular, if the 
integral in the y-direction of the azimuthal component of the 
EMF at each radial face is not equal, then the net vertical flux 
in the domain is not conserved. Since the dynamics of the MRI 
is sensitive to the amplitude of the net vertical flux, this can rep- 
resent a serious problem. To circumvent this issue, we remap 
the azimuthal component of the emf at each radial face, and use 
the arithmetic average of computed for each grid cell at each 
radial face, and the remapped value of 8y from the correspond- 
ing grid cell on the opposite face, to update the magnetic field 
in the cells next to the boundary. Since the remap conserves the 
integral of in the y-direction, the integral of the averaged 8y 
used to update the field on each boundary will be identical, and 
so this procedure conserves the net flux in the vertical direction 
to machine precision. In principle, this same procedure could 
be used for the fluxes of all the conserved quantities, to restore 
conservation in the domain. In particular, it can also be used 
to remap so that the volume averaged azimuthal flux of the 
magnetic field is conserved to machine precision. In practice. 



we find that with the orbital advection algorithm described in 
§5, the hydrodynamic variables are conserved accurately even 
without this extra step, so that we generally only correct in 
this fashion. 

In the same way that the net vertical flux in the domain is not 
conserved without a special treatment of the remap of E^, the 
net toroidal flux is not conserved without a special treatment of 
the remap of We have run a variety of simulations of the 
MRI in three dimensions, and find the net toroidal flux intro- 
duced in this way is very small: the magnetic energy associated 
with this flux shows random fluctuations with peak amplitude 
of around 10"^ of the thermal pressure for a zero-net- vertical 
flux simulation with an initial /3 = 4000. However, we have 
found that matching the order of the reconstruction used in the 
ghost cells with the order used in the integrator, that is using 
third-order reconstruction in the ghost cells when third-order 
reconstruction is used throughout the rest of the domain, is im- 
portant to keep this anomalous flux small. Using second-order 
reconstruction in the ghost zones produced a mean toroidal flux 
that is still very small, but 10 times larger than when third-order 
reconstruction is used in the ghost zones. 

Efficient use of a large number (> lO-*) of processors on mod- 
ern parallel computers requires domain decomposition of the 
grid in all three coordinate directions. However, paralleliza- 
tion of the shearing box boundary conditions can require com- 
plex book-keeping if domain decomposition in the azimuthal 
direction is allowed, because the data in the ghost zones for 
each processor on the radial boundaries can require communi- 
cation with up to 3 processors on the opposite boundary, de- 
pending on the configuration of the domains as they shear. To 
simplify the implementation, we first apply periodic boundary 
conditions (which are easily parallelized) to pass ghost cells be- 
tween processors in the radial direction, followed by a remap in 
the azimuthal direction. The latter is easily parallelized using a 
cyclic shift operation across those processors that store neigh- 
boring grid patches in the azimuthal direction. The cyclic shift 
is also easily parallelized, and always occurs between the same 
set of processors along the boundary. The disadvantage of this 
method is that it adds extra MPI calls, however this is more than 
offset by the reduction in complexity of the code, and the fact 
that the resulting communication pattern is more regular and 
therefore in most circumstances likely to be more efficient. 

5. ORBITAL ADVECTION 

Orbital advection methods (Masset 2000; JGG; Johansen et 
al. 2009) exploit the fact that the background orbital motion v^- 
is time-independent and varies only in radius (the it-coordinate 
in the shearing box). Thus, by decomposing the velocity into 
two parts 

v = Vjf + v' (41) 
where the components of the velocity fluctuation vector v' are 







Vx 






Vy + q^QX 
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the equations of MHD in the shearing box (equations[T|through 
lU can be rewritten as 



dp dp ^ ^ ,^ 



:0, 



system equations |53] through |56l are just the usual equations of 
MHD, but written in terms of a vector of conserved variables 
q = (p,pv',£'',B), the standard CTU+CT algorithm in Athena 
(4^n be used to integrate these equations, with the appropriate 
Crank-Nicholson time differencing of the shearing box source 



^^ + Vjf^^ + V- [pv'v' + T] =21]o(pv^)i + (g-2)rio(pVt)j-P^^™^ ^" momentum fluctuation equation, and the gradi 
dt cly gjjl; of the effective gravitational potential in the energy equa 



- + VK^ + W- [E'y' + T . v'] = py' • V<i>^3 + - 



dt 



I /, 9vjf lion. These source terms are in fact simpler than in the origi 



SB 

'dt' 



■ V X (vjf X B) - V X (v' X B) = 0, 



where the total energy E' contains the kinetic energy in the ve- 
locity fluctuations (but not the background shear flow) 
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(47) 



and the effective potential for the shearing box in the orbiting 
frame contains only the vertical gravity 

1 



2„2 



SB ■ 



(48) 



Note the second term on the LHS in each of equations [43] 
through |46] is a linear advection term with characteristic speed 
x'K- This suggests a numerical algorithm based on splitting these 
equations into two systems, the first containing only the linear 
advection operators 



dt dy 



■0, 



dpy' dpy' 

dE' 
Ik 



dy 
(9B 

- — Vx(v/rxB): 
ot 



dy 
dE' 

+ VK-^ = {B^By-pv'yy) 



dVK 

dx 



0. 



(49) 



(50) 



(51) 



(52) 



Note the extra terms on the RHS in the energy equation. These 
terms cannot be written as flux gradients, and must be treated 
as source terms. We discuss them in more detail below. 

The remaining system of equations are the usual equations 
of MHD, written in terms of the velocity fluctuations v' rather 
than the velocity v, that is 



0, 



(53) 



dpy' 
"dT 



+ V- [py'y' + T]=2no(pv[:)i+(q-2)no(pv'Sj-pfll0i) 



dE' 

~dr 



+ V- [£V + T-v'] =pv'-V$s 



SB J 



f-Vx(v'.B) 



0, 



(55) 



(56) 



Note that since these equations written using pvj, and E' as 
the conserved variables, instead of p\\. and E, the shearing 
box source terms are modified compared to the original system 
equations [T] through |4] 

Developing numerical algorithms to solve these two systems 
of equations is straightforward. In particular, since the second 



system, and the techniques described in §3 apply directly. 
Since the variables being updated are the momentum fluctua- 
(4^Jons themselves, the source terms added to the reconstruction 
step and transverse flux gradient corrections (see §3.2) are mod- 
ified appropriately. Moreover, since the Riemann solver now 
returns the fluxes of the momentum fluctuations directly, the 
conversions in equations [3T] and [32] are no longer necessary. 
Note that the CFL stability condition for this system is based 
only on the amplitude of the velocity fluctuations, that is 

5x 5y 5z 



St = Co min 



( K I + Cf^u,^ ' ( I y; I + Cy,)/,M ' ( I I + Cf^ijM 



(57) 

where Co < 1 is the CFL number, and C/a, C/x, C/.i are the fast 
magnetosonic waves speeds in the x-, y-, and z-directions re- 
spectively, and the minimum is taken over all grid cells. Since 
the background orbital flow becomes supersonic for \x\ > H/q, 
where H = CJHq is the scale height in the disk, the time step 
can be much larger using orbital advection in domains which 
span roughly H or larger in the radial direction. This is the 
primary advantage of the method. 

Since the first system of equations (|49]throughl52Ti are linear 
advection operators in one dimension, the numerical algorithms 
to integrate these equations are particularly simple. The finite 
volume discretization of the first three of these equations can be 
written as 



U 



n+l 



■u" 



'i,j,k - "^iJM-i^f /^y)(fi.j+l/2,k- fi,j-l/2,k) 

where U represents each of the conserved quantities, and 
fij+i/ik the upwind fluxes of these quantities at cell interfaces 
in the y-direction. These fluxes are simply the integral of U 
over the domain VKt upstream of the appropriate interface, for 
example 



(58) 



U''{xi,y,Zk)dy 



(59) 



(Vj-i/2)-i'*;i5; 



Numerically, this integral is converted into a finite sum over 
all the grid cells upstream of the interface in the y-direction. In 
general, the domain of dependence VK^it will span a non-integer 
number of cells. For those cells entirely contained in the do- 
main, the volume-averaged value of U is added to the sum. For 
the one cell which is only partially contained within the inte- 
gral, then the higher-order reconstruction functions described 
in §4 of SGTHS are used to evaluate the fraction of the con- 
served variable that contributes to the integral. Note there is no 
time step constraint for stability of this algorithm. The domain 
of dependence can be arbitrarily large, and span any number of 
grid cells. The step simply represents a conservative remap of 
the solution in the j-direction by a distance VK^t. 
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There are two aspects to the orbital integration algorithm that 
warrant some discussion. The first is the integration of the in- 
duction equation. Note we have written equation |52] in a form 
that suggests the use of the CT algorithm with an effective emf 
given by £ = -\k x B. The components of this emf are 

{E,,Ey,£d = (-VKB,,Q,VKB,) (60) 

The centering of the components of the emf used in Athena 
are shown in figure 1 in SGTHS. The discrete form of the CT 
update for each component of the magnetic field is given by 
equations 16 through 18 in SGTHS. The CT algorithm for or- 
bital advection simply requires the calculation of the effective 
emf by integration of each component of the electric field £ 
upstream of the appropriate cell edges distance vkSi, that is 

n+l/2 p'i-^n- 

^".i.j-\/Xk-il2 = ~ VKB".(xi,y,Zk-i/2)dy (61) 

n+l/2 

Ki-i/2.j-i/2,k= VKB'^{xi.i/2,y,Zic)dy (62) 

tyj-i/2)-''*;<5' 

By using a CT discretization of equation |52] to evolve the 
magnetic field in the orbital advection step, we preserve the 
divergence-free constraint to machine round-off during the or- 
bital advection step. Moreover, note that there are no source 
terms required in equation |52] The growth or decay of By due 
to the shear is naturally captured in the CT difference formu- 
lae, avoiding any divergence-free interpolation as required in 
the algorithm of JGG. 

The second aspect of the orbital advection step that requires 
further discussion is the integration of the energy equation ISTI 
Note the source terms that appear on the RHS, which represent 
the work done by Reynolds and magnetic stress due to the radial 
shear of the orbital motion. The finite-volume discretization of 
equation |5T] requires a time- and volume-averaged approxima- 
tion for the source terms. In the Lagrangian frame (comoving 
with the fluid during the remap), all quantities except By are 
constant, and the time evolution of By in the Lagrangian frame 
is particularly simple, By(t + 6t) = By{t) + By(dvK/dx)6t, where 
dvK/dx = -qVl() This suggest a "Lagrangian-step-plus-remap" 
algorithm for the energy equation. In the Lagrangian step, the 
total energy is updated using the time-averaged source terms 

Kjjc = E;jk-qnoStB,jj,k(Byjj40)-B,^ij,kqnQSt/2) 
+ ^rjo^f K,,j-,i^,'v.>,* ^^^^ 

where B^jj^k and Bv,!j,/t(0) are the volume-averaged, cell- 
centered components of the magnetic field at the start of the 
orbital advection step. The remap step then uses equation |58] 
applied to the total energy with the source terms added, j i^, 
to complete the update of the total energy. Comparisons of cal- 
culations of the MRl with an adiabatic equation of state, both 
with and without orbital advection (see next section), show that 
this algorithm for the energy source terms works well. Note 
that the only way in which the total volume averaged energy in 
the domain can change is via a non-zero total stress at the radial 
boundries. Additional tests have shown our method conserves 
the total energy if the stress at the boundaries is zero. 



5.1. Tests of Orbital Advection 

To test our orbital advection algorithm, we have run calcu- 
lations in both two dimensions (in the x-y plane), and full 
three dimensions. Our calculations span < x < L^/l, 

-Ly/2 <x< Ly/2, and -Lj2 <z< Lj2. We set the orbital fre- 
quency = 10"^, the shear parameter q = l>/2, and unless oth- 
erwise stated adopt an isothermal equation of state with sound 
speed Cj (we give the value of C, used for each test below). 
As before, we set the initial density po = 1 and orbital velocity 
Vy.o = -q^QX, and we use third-order reconstruction, the HLLC 
(for hydrodynamics) or HLLD (for MHD) Riemann solver, and 
the CTUh-CT unsplit integrator for all the tests. 

Our first test of orbital advection is the evolution of a hy- 
drodynamic shearing wave (Johnson & Gammie 2005, Balbus 
& Hawley 2006; Shen et al. 2006). We use a domain of size 
Li = Ly = 1, an isothermal equation of state with sound speed 
Cj = 1.29 X 10""^, and an initial wavevector 27rko/L = (-8,2). 
We expect the maximum amplification of the wave amplitude 
to be 17, at time ffio = 8/3. Figure 5 compares the amplitude 
of the kinetic energy associated with the x-component of the 
velocity both with and without the orbital advection algorithm 
at various resolutions. In both cases, the evolution is accurately 
captured provided there are at least 8 grid points per wavelength 
in the initial conditions. Most importantly, there is virtually no 
aliasing in the solutions at late times. In both cases, the kinetic 
energy decays over 5 orders of magnitude, with no periods of 
sustained growth after the maximum. 

Our first full MHD test of orbital advection is the advection 
of a weak (/3 = 2 x 10^) field loop. This test has been useful for 
developing the CT algorithm implemented in Athena (Gardiner 
& Stone 2005a; 2008; Stone et al. 2008). We use a domain 
of size Lx = 3 and Ly = 8, an isothermal equation of state with 
Cj = 10"^, and initialize the field loop following the method de- 
scribed in GS05 centered at the origin. To make the problem 
more complex, we introduce a uniform radial velocity Vv = C,, 
which causes the loop to execute epicyclic oscillation as it is 
sheared. The amplitude of these oscillations is large enough so 
that the outer one quarter of the loop crosses the radial boundary 
at the extrema of the motion. This tests whether our implemen- 
tation of the shearing box boundary conditions described in §4 
can maintain the integrity of the loop. Figure 6 shows plots of 
the field lines at four times in the evolution. Note there are no 
indications of the stripes which appear if the EMFs are com- 
puted incorrectly (see figures 2 and 3 in GS05), confirming the 
upwind CT method developed by GS05 to compute the emfs 
also works well with orbital advection. In addition, there are no 
features in the loop associated with the boundaries, indicating 
our implementation of the shearing box boundary conditions is 
accurate. 

Another sensitive MHD test recently introduced by JGG is 
the evolution of a compressible shearing wave (Johnson 2007). 
We have repeated the test shown in figure 1 1 of JGG, and com- 
pared the resulting solution to a numerical integration of the 
ODEs that describe the analytic solution to the problem kindly 
provided to us by B. Johnson. The calculation uses a domain 
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of size Lx = Ly = = L = 0.5, and a numerical resolution of 
A' X N/2 X N/2, where N = 16, 32, and 64 (these values are 
chosen to match those used by JGG). An isothermal equation 
of state with C, = 1 is adopted, and for this test we use llo = 1 . 
Although the details of how to initiahze such tests are given in 
the Addendum to JGG, to be specific we give the initial con- 
ditions we actually used in Athena for this test, namely den- 
sity p = 1 + 5 .48082 X 1 0"'' cos (ko • x), velocity fluctuations v' = 
(-4.5856 X 10-^ 2.29279 x 10-^ 2.29279 x 10-*)cos(ko -x), 
and magnetic field B/V4^ = (5.48082 x 10-^1.0962 x 
10-^,0)cos(ko •x) + (0.1,0.2,0), where the wavevector ko has 
components (-2, 1, l)(27r/L). Figure 7 shows the time evolu- 
tion of the volume-averaged perturbations in the real part of the 
azimuthal component of the magnetic field, defined as 

6By = [ 2(B, - B,) cos (k(f ) • x) (64) 
Jv 

where By is the numerical solution computed by Athena, By = 
0.2-0.150o^ is the expected evolution of the real part of 
the field, and the wavevector has components k(f ) = (-2 + 
q^lot,!, \)(2tt/L). Note the highest resolution solution, which 
has 32 grid points per wavelength in each direction initially, is 
indistinguishable from the semi-analytic solution. 

It is instructive to compare the nonlinear saturation of 
the MRl in calculations both with and without orbital ad- 
vection. We present the evolution of two different three- 
dimensional calculations, one with orbital advection and the 
other without. Both start with no net vertical flux, B = 
(0,0,(2po//3)'/^sin27rx/LJ, where (3 = 400, use a domain of 
size Lx = Ly = SH, and L^ = H (where H = Cs/ f2o = 1), and adopt 
an isothermal equation of state with C, = 10"^. The resolution 
is 32 /H in each dimension. Figure 8 plots the time evolution 
of the magnetic energy in both cases. In the nonlinear regime, 
the timestep in the calculation with orbital advection was on 
average 4.3 times larger than the calculation without. The re- 
sults show the solutions are in excellent qualitative agreement. 
Because the nonlinear evolution of the MRI is chaotic (Winters 
et al. 2003), we do not expect exact agreement, but only the 
same values for suitably time-averaged quantities. The level 
of agreement between our two calculations is similar to that 
reported by JGG. These authors point out that because the trun- 
cation error varies with x in the shearing box without orbital 
advection, there are unphysical patterns in the radial profile of 
time-averaged quantities such as the stress and density. Fig- 
ure 9 plots these profiles for this calculation. Note that without 
orbital advection, we observe a minimum in the density, and 
a maximum in the stress, near the center of the calculation. 
This is likely caused by a minimum in the truncation error in 
the region of this point, since the velocity amplitude (including 
the radial, azimuthal, and vertical components) is a minimum 
near there. However, with orbital advection the profiles of the 
density and stress are clearly much smoother; they only show 
small ampUtude variations associated with the MRI-driven tur- 
bulence. 

These calculations do not use the corrections to the shearing 
box boundary conditions (Gressel & Ziegler 2007) that are re- 
quired to conserve mass and momentum exactly (energy is not 



conserved with an isothermal equation of state). Nonetheless, 
we find that at the end of the calculation without orbital ad- 
vection, mass is conserved to one part in 10^, and with orbital 
advection it is conserved to one part in 10^. Since these simu- 
lations are only run for 16 orbits, we have also checked mass 
conservation at various resolutions for zero net vertical simula- 
tions in a domain of size = H, L,. = 4H, and = H, using 
resolutions of 32/H, 64/H, m/H, and 256/i/. This study 
was conducted to confirm the lack of convergence of the am- 
plitude of the Maxwell stress in zero-net flux unstratified boxes 
discussed by Fromang et al. (2007), and indeed we confirm this 
result. A variety of other results from this study will be reported 
elsewhere. Here we note that without orbital advection, mass 
conservation can be a problem at the lowest resolutions. At 
32/H, roughly 8% of the initial mass is lost by 100 orbits, with 
the mass loss rate essentially constant with time. The mass loss 
rate decreases with resolution, however even at 256///, 1.5% 
of the mass is lost by 100 orbits. The use of orbital advection 
greatly improves conservation. After 100 orbits, we find only 
0.03% of the initial mass is lost by 100 orbits with orbital ad- 
vection, and this number is independent of resolution. 

Finally, we have also used simulations of the MRI with an 
adiabatic equation of state to investigate the implementation 
of the source terms in the energy equation in the orbital ad- 
vection step. We have repeated a simulations with net ver- 
tical flux B = (0,0,(2po/(3y/^), where /3 = 400, in a domain 
of size Ljf = H, Ly = 2ttH, and = H, using a resolution of 
64 X 128 X 64. This is a repeat of run Z2 in Gardiner & 
Stone (2005b). The calculation starts with po = 10~^, and uses 
7 = 5/3. We have run calculations both with and without orbital 
advection, and compared the resulting time evolution of the to- 
tal energy with figure 3 in Gardiner & Stone (2005b). Again, 
because of the chaotic nature of the MRI, we do not recover 
exact agreement in the two calculations, however both show 
growth of the total energy at very similar time-averaged rates, 
with aU of the increase after saturation occurring in the internal 
energy density. 

6. APPLICATION TO THE MRI 

The test results presented in the previous section demonstrate 

the orbital advection algorithm is at least as accurate (although 
perhaps not any more accurate) than integrations that do not 
use it. The primary advantage of orbital advection, however, is 
not that it is more accurate, but that it removes the background 
shear flow from the time step stability limit, and therefore en- 
ables much more efficient studies of accretion flows over a wide 
range of radii. Moreover, the new method based on CT intro- 
duced in this paper is simpler than previous approaches, and 
preserves the divergence free constraint to machine precision. 

Figure 10 shows images of the density and azimuthal velocity 
fluctuations in a shearing box simulation of the MRI in a very 
wide domain of size = Ly = 32H, and = H, using a resolu- 
tion of 32/H in each dimension, or 1024 x 1024 x 32. The ini- 
tial conditions and parameters of the calculation are identical to 
that presented in figure 8. Interestingly the spiral density waves 
excited by the MRI are immediately obvious, and have the same 
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pitch angle as in smaller domains, indicating the box size does 
not determine these features (Heinemann & Papaloizou 2008a; 
b). Recently Johansen et al. (2008) have reported long-lived 
density features in domains with wide vertical extent, and a 
more careful analysis of this and other runs reproduces this re- 
sult. The time averaged Maxwell stress in this calculation is 
little different than the value in much smaller domains, indicat- 
ing the turbulent stress must be localized on scales ^ H (Guan 
et al. 2009). 

JGG found unphysical features in the time-averaged radial 
density profile even with orbital advection, wherever the orbital 
shear displacement was close to an integer number of zones. 
In figure 11 we plot this profile for the calculation shown in 
figure 10. Using the time step measured from this simulation, 
the radial locations where Stq^o/5y is an integer are plotted as 
vertical dashed fines. There are no significant extrema at these 
locations, indicating the stress and truncation error are smooth 
with radius. This is another advantage of the MHD orbital ad- 
vection algorithm developed in this paper. 

In summary, we have described the inclusion of source terms 
for the shearing box approximation in the Athena MHD code. 



including a Crank-Nicholson time differencing that preserves 
the ampUtude of epicyclic oscillations exactly. We have also de- 
scribed an orbital advection algorithm based on CT for evolving 
the induction equation to preserve the divergence free constraint 
on the magnetic field. We have shown this algorithm provides 
more accurate solutions at less computational cost. These al- 
gorithms have already been used to study hydrodynamic turbu- 
lence in the shearing box (Shen et al. 2006). A number of new 
studies of the MRI in wide radial domains in both unstratified, 
and vertically stratified disks are underway (Davis et al. 2010). 

We thank Sebastien Fromang, Charles Gammie, John Haw- 
ley, Bryan Johnson, E. Ostriker, and Jake Simon for helpful 
discussions, and S. Fromang for providing the initial condi- 
tions for the nonUnear density wave test, and B. Johnson for 
providing the reference solution for the MHD shearing wave 
test. Financial support was provided by the DOE through DE- 
FG52-06NA26217. Simulations were performed on computa- 
tional faciUties at the Princeton Institute for Computational Sci- 
ence and Engineering, and through resources provided by NSF 
grant AST-0722479. 
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Fig. 1 . — Time evolution of radial velocity amplitude, normalized by the initial value, in epicyclic motion in a box of radial extent Lx-\ (solid line), = 10 
(dashed line), and Ly: = 50 (dotted line) computed on a grid of 64^. There is never amplitude error in the solutions, and only small phase error in the largest box (in 
which the timestep is about 0. 1 orbits). 
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Fig. 2. — (Leji). Fractional change in density in an axisymmetric nonlinear density wave at f = lOr (dotted line), 50r (dashed line), and lOOT (sohd line), where 
T = 5254.31 is the wave period, on a grid of 40 cells. (Right). Fractional change in the wave amplitude as a function of time. Both quantities are defined in the text. 
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Fig. 4. — Time evolution of 0.5BI in tlie axisymmetric MRI. The solid lines are for third-order spatial interpolation, the dashed lines are for second-order The 
numerical resolution of each pair of lines is labeled. 
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Fig. 5. — Time evolution of tlie radial velocity amplitude in an incompressible shearing wave at numerical resolutions of 32^, 64^^, and 128^, corresponding to 4, 
8, and 16 grid points per wavelength initially. Solutions computed both without orbital advection (left), and with orbital advection (right) are shown. The dashed 
line in both cases is the analytic solution, and for each curve the amplitude is normalized by the initial value. Particularly notable is the extremely small level of 
aliasing. 




Fig. 6. — Field hues of an initially circular field loop advected by epicychc motion and sheared by the background flow, computed using orbital advection. Twenty 
field lines are shown at each times of (from left to right) t = 0, 0.3, 0.7, and 1 orbit. 
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Fig. 7. — Time evolution of the real part of the wave amplitude in the azimuthal component of the magnetic field for an MHD shearing wave test, computed using 
orbital advection. The solid Unes are at a resolution of 8, 16, and 32 grid points per L^, while the dashed line is the semi-analytic solution. 
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Fig. 9. — Radial profile of the average density (left) and stress (right) in a zero-net flux MRI calculation computed with (dashed line) and without (solid line) 
orbital advection. The averages are taken over all y and z, and over orbits 7-16 for the solid line, and 70-100 for the dashed Une. 
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Fig. 10. — Images of the density (left) and azimuthal velocity fluctuations scaled to the sound speed (right) from a zero-net flux MRI calculation in a box with 
radial dimension L,- = 32H at ; = 1 6 orbits, computed with orbital advection. 
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Fig. 1 1 . — Radial profile of the average density in a zero-net flux MRI calculation in a box with radial dimension Lv = 32H. The average is taken over all )' and z, 
and over orbits 7-16. For the typical timestep used in the calculation the azimuthal displacement in the orbital advection algorithm is an integer number of zones at 
the locations marked by the vertical dashed lines. There is no evidence of density minima at these locations. 



